Shock waves in the dissipative Toda lattice 
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^ . Abstract 

. ^ We consider the propagation of a shock wave (SW) in the damped Toda lattice. 

X ■ The SW is a moving boundary between two semi-infinite lattice domains with dif- 

ferent densities. A steadily moving SW may exist if the damping in the lattice is 
represented by an "inner" friction, which is a discrete analog of the second viscosity 
in hydrodynamics. The problem can be considered analytically in the continuum 
approximation, and the analysis produces an explicit relation between the SWs 
velocity and the densities of the two phases. Numerical simulations of the lattice 
equations of motion demonstrate that a stable SW establishes if the initial velocity 
is directed towards the less dense phase; in the opposite case, the wave gradually 
spreads out. The numerically found equilibrium velocity of the SW turns out to be 
in a very good agreement with the analytical formula even in a strongly discrete 
case. If the initial velocity is essentially different from the one determined by the 
densities (but has the correct sign), the velocity does not significantly alter, but 
instead the SW adjusts itself to the given velocity by sending another SW in the 
opposite direction. 



1 Introduction 



The analysis of dynamical behavior in linear and nonlinear lattices is of obvious interest 
for various branches of solid state physics ]I[ @| . It is known that anharmonic interactions 
in the lattice may produce stable localized collective excitations in the form of solitons 
IfJ, [|] in the absence of dissipation. The dissipation naturally damps solitons M, although 
they may be supported, in certain cases, by an external AC drive ||. 

The aim of the present work is to consider a different nonlinear wave in the damped 
lattice, namely a shock wave (SW) or kink that is realized as a localized transient region 
between two parts of the lattice with different spacings, i.e., different densities. Evidently, 
the shock may be produced by an initial deformation of the lattice. Another way to 
produce a shock is to apply a strong short pulse to an edge of the lattice (e.g., by an 
impact of another body). Mathematically, the latter type of the shock corresponds to 
the von Neumann problem: in an undeformed lattice, the particles with the positive and 
negative numbers are given initial velocities ± Vq @ • A similar but different problem was 
considered by Kaup for the semi-infinite lattice, in which the first particle is driven at 
a constant velocity, while others are initially at rest. 

For the non-dissipative (exactly integrable) Toda lattice (TL) model, evolution of the 
von Neumann shock has been analyzed in a rather full detail ||. It has been demonstrated 
analytically and numerically that the evolution is essentially non-stationary, producing a 
lot of oscillations behind the shock. (This seems to be qualitatively similar to the decay 
of an initial step configuration in the Korteweg - de Vries (KdV) equation, which is a 
continuum limit of the TL ||.) The principal difference between the von Neumann's 
problem and the problem to be considered in this work is that the former one is initially 
dealing with a homogeneous lattice and infinitely large kinetic energy (if the system is 
infinite), while we want to consider a lattice that is initially deformed so that it has an 
infinite reservoir of potential energy. 

The presence of suitable dissipation in the lattice seems to change the behavior in an 
essential way: With dissipation it may be possible to obtain a steadily propagating shock 
between the phases (lattice domains) with different densities. It is interesting that the 
inclusion of friction, which renders the problem physically relevant, also is necessary for 
the stable propagation of certain waves. 

In this work we will consider shock propagation in the framework of the damped TL 
although the general results obtained below should remain qualitatively correct for a fairly 
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broad class of nonlinear dynamical models. In the TL one can introduce two different 
types of friction, which can be naturally called outer and inner. The model with the outer 
friction is based on the following lattice equation of motion: 

where IS cl> displacement of the rath particle from its equilibrium position. The dissi- 
pative term in Eq. (1) implies friction between the lattice and some external substrate, 
hence the name "outer" friction. However, for applications to solid state physics, it seems 
more natural to consider friction of purely internal origin, which is similar to the second 
viscosity in hydrodynamics, i.e., resistance against change of the density of the lattice. 
The TL model with inner friction is based on the following equation of motion: 

x n + 7 (2x n - x n _! - x n+1 ) = e -(*»-*»-0 - e -(*»+i-»n) . ( 2 ) 

A local density of the lattice is determined by the local relative displacement r n = x n — 
x n _i, which in the inner friction case obeys the equation 

r n + 7 (2r n - r„_i - r n+1 ) = 2e~ r " - e' Tn+1 - e^- 1 . (3) 

The variables x n in Eqs. (1) and are not necessarily full coordinates of the particles, 
but may be, as it was mentioned above, displacements of the particles from some regularly 
spaced (with distance a) equilibrium positions. The full distance between the particles 
will then be a + r n , and the local density of the lattice (a + r n ) . This also implies that, 
in what follows below, r n may take both positive and negative values. 

In this paper our main concern is the study of SW's governed by Eq. (|3]). In the 
second section we obtain some analytical results using the continuum approximation. We 
can demonstrate, e.g, that a steadily traveling SW always exists, and although it cannot 
be generally found in an explicit form, we obtain an analytical expression for its velocity, 
which turns out to depend on the asymptotic values of r in front of and behind the SW, 
but not on the friction coefficient 7 in Eq. (|3|). This velocity is sandwiched between the 
sound velocities of the two phases, and the motion is towards the less dense phase. In 
the third section, we display results of numerical simulations of Eq. ([3]). When the initial 
velocity of the SW has the correct sign, we always end up with a stably propagating wave 
of a permanent shape, but when the direction is wrong the SW gradually flattens out. 
We find also that when the SW is started with a wrong magnitude of the velocity (but 
with the correct sign), the SW readjusts its height so that a fairly good agreement is 
reached with the analytically predicted formula. The change of the height is carried out 
by emitting an additional SW in the opposite direction 
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2 Analytical consideration 



As the first step in analyzing the possibility of a steady propagation of the SW, it is 
natural to define the "quasi-momentum" 

ra=+oo 

p= £ ( 4 ) 

n=—oo 

Summing up all the equations (1) over n from — oo to +oo, one immediately obtains the 
equation 

f + 7 P = 0. (5) 

Since the friction coefficient 7 is positive, Eq. (f|) tells us that P — > as t — > 00. On the 
other hand, it is obvious that a steadily propagating SW must have a nonzero value of 
the quasi-momentum. Therefore, a steadily propagating SW is not possible in the model 

(!)• 

Proceeding to the model (|]) with the inner friction, one notices, first of all, that the 
same trick yields, instead of Eq. the equation ^ = 0. Thus, in this case a constant 
nonzero value of the quasi-momentum is possible. 

Let us now consider the continuum approximation. Assume, as usual, that r n in (Q) 
changes appreciably only on a scale much larger than the lattice spacing. One may treat 
n as a quasi-continuous variable (coordinate), and then in the lowest approximation, Eq. 
(H) goes over into the continuum equation 

r tt - ir tnn = - (e' r ) , (6) 

where the subscripts stand for the corresponding partial derivatives. Note that, generally 
speaking, we do not assume the value of r n small, and therefore we do not expand the 
exponential. For the time being we do not take into account the fourth derivative on the 
right-hand side. In what follows below it will be added to Eq. (|||) when necessary. 

One should now look for traveling- wave solutions to Eq. (||) in the form r = r(z), z = 
n — Vt, V being the wave's velocity. Substituting this into Eq. (J6|), one arrives at an 
ordinary differential equation which can be integrated twice. With regard to the condition 
that we are interested in solutions which remain finite at n — > ±00, the final form of the 
integrated equation is 

lV J = - (V 2 r + e~ r - C) (7) 

C being an arbitrary constant of integration. In what follows, we will consider the asymp- 
totic values r± = r(z = ±00) as given parameters (as it was explained above, they are 
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determined by the density of the lattice in front of and behind the SW), and since the 

dr 

dz 



derivative must vanish at the two infinities, we find 



V 2 = (r+ -r^Y 1 (e- r ~ - e~ r+ , . 



C = (r+ - r^y 1 (r + e~ r - - r_ e - r +) . (9) 

Note that the friction coefficient does not enter in these relations. For the maximum value 
of the slope we get 

max (9 = ^ 2iogi/2 - i/2+c) ' <io) 

this value is reached when r = — log\^ 2 . 

The linearization of Eq. (|5|), with 7 = 0, around an arbitrary constant value r (r = 
tq + p, p being a small variable part of r) yields a D'Alembert equation for p, which 
determines the sound velocity s in the homogeneous lattice with r = r : 

s 2 = e~ r °. (11) 

It is easy to check that the shock's velocity given by Eq. (|8|) always lies between the 
sound velocities corresponding to r — r±, and in the limit r + — r_ — > coincides with 
the corresponding sound velocity. Furthermore, the velocity (§) coincides with the sound 
velocity at the maximum slope (|T0|). 

An explicit profile of the SW cannot be obtained from Eq. ([?]) in a closed form. 
However, this can be done in the limiting case |r + — r_| r±. In this case, Eq. (^) 
reduces to 

dv 

^is- 1 — = (r+ -r)(r-r_), (12) 

where s is the sound velocity ([11]) corresponding to r± (in this approximation, we neglect 
the difference between them). The solution to Eq. (jl2|) is a typical kink 



r(z) = — [(r+ — r_) + (r+ — r_) tanh^] , (13) 

where C = |s7 _1 (r + — r_)^. 

Let us now briefly consider the applicability conditions for the approximation consid- 
ered. As it follows from the explicit solution (|12|) , the size L of the SW can be estimated 
as follows: 

L ~ -y/V5r (14) 

where 5r = r + — r_. Thus, the basic condition for application of the continuum approxi- 
mation, L ^> 1, takes the form 

7 > V5r. (15) 
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Above, we have neglected the fourth derivative on the right-hand side part of Eq. 
(0). This is justified if, expanding the exponential, one has the second derivative of the 
squared variable part p of r much larger than the fourth derivative p nn nn- Using again 
the solution fll^D, one concludes that the corresponding condition is 

7 2 > V 2 6r. (16) 



In the general case, when 5r is not specially small, the conditions (|15|) and (|T6| ) are actually 
equivalent, i.e., there is no necessity to add the fourth derivative to Eq. (|6|). However, in 
the case of the weak SW, when Sr is small, the inequality (|I5|) is less restrictive than (|16"D. 

Representing, as above, r in the form r + p, one obtains an equation for the (small) 
variable part p. This equation can be immediately integrated once, so that one ends up 



with the known Korteweg-de Vries - Burgers equation | 10[| : 



1 . . 

2p t - Ipzz ~ spp z + —Spzzz = , (17) 

where z is the same traveling coordinate as in Eq. (|7|), and s is the sound velocity (|TI|) cor- 



responding to r = ro- In the work [T(| it has been demonstrated that Eq. (|17D has steady 
traveling- wave solutions in the form of a SW with an oscillating tail (it was approximately 
described in [9] as a bound state of several KdV solitons). This solution represents the 
lattice SW in this limiting case. 

Also in this case an exact analytical solution for the shock's shape is not available, 
but one can easily find its velocity from Eq. (|r7j). Going back to the original coordinate 
n, one can eventually obtain the full shock's velocity as follows: V = s — \s(p + + p_), 
where p± = r± — r . On the other hand, expansion of the general formula ([§D for small 
r + — r_ yields exactly the same result. Thus, the addition of the fourth derivative affects 
the shape of the SW, but not its velocity. 

To conclude this section, it is relevant to note that the continuum approximation has 
another limit of applicability when the lattice density in front of the SW becomes too 
small, i.e., r + is very large. Formally, Eq. (0) predicts that the velocity vanishes in this 
limit. However, it follows from Eq. (|7]) that the size of the corresponding SW is, in the 
continuum approximation, L ~ r yVr_e r ~. Thus, L shrinks when V vanishes, and the 
continuum approximation is no longer valid. 
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3 Numerical results 



We have simulated Eq. ([|) numerically in the logarithmic form 

d 2 ln(l + v n ) + -ydt [2 ln(l + v n ) - ln(l + v n+1 ) - ln(l + v n ^)} = v n+1 + v n _ x - 2v n , (18) 

where v n = exp(— r n ) — 1, or more precisely, in its two-component form 

d t v n = (z n - i n+ i)(l + v n ), (19) 
d t i n = v n -i -v n - 7(2i n - i n+1 - z n _i)- (20) 

System (19 j20|) is more suitable for numerical integration than the original one, because 
there is no need to evaluate transcendental functions, and therefore local numerical errors 
can be significantly reduced. We have used the Bulirsch-Stoer algorithm as the integration 
method: the total number of the lattice sites was 2000. With a lattice of this length it 
was possible to study long-term evolution of the initial state. 

As the initial state we took f n (0) where v n (t) is a SW of the form 

Vn (t) = v Q + ^Av (1 - tanh[A;(n - n - Vt)}) , (21) 

and then i n (0) was obtained from (^H) with (pT|) and 7 = by integrating in time: 



i (0) = l n ( cosh[fc(n-l-n )] \ 

UUj 2Vk ln { cosh[k(n - n )} | ^ 

[Note that i n is only needed up to a constant additive part.] In these formulae the velocity 
V and slope k are free, but it turns out, that if we want the best initial value we should 
take V from (||) with r + = — ln(l + v + At>) and r_ = — ln(l + v ), and the slope of the 
initial state from (|Hj) (recall r n = — log(v n + 1)): 

k = --^(vHnV 2 -V 2 + C) . (23) 
At>7 v ' 

The time evolution of a typical solution is shown in Figs. 1(a) and 1(b) in the form of 
snapshots of the lattice at various times. If the dissipation is zero (Fig. 1(a)) the front of 
the SW starts immediately to produce oscillatory ripple waves at its top shoulder. These 
oscillations increase in amplitude and extend slowly to the left of the shoulder. If we 
have, however, moderate dissipation (7 = 1 in Fig. 1(b)), the SW travel with constant 
shape and velocity without any oscillating tail. We have varied the amplitude parameters 
Vq and Av from 0.5 to 10, and the dissipation factor 7 from 0.1 to 2 (which implies that 
the width of the slope varies from 2 to 50 lattice points) and always obtained similar 
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results. Typical final SW profiles are presented in Fig. 2 using three different values for 
the dissipation factor. The velocity of the SW can be predicted extremely well (deviation 
is less than 0.2%) using the continuum limit result (|j), even when the system is highly 
discrete. 

The convergence of the initial solution towards the final steady state is illustrated in 
Fig. 3(a). This figure contains several snapshots of the lattice plotted on top of each other, 
after a certain shift relative to the lattice. To determine the shift, we first sought for the 
position of the SW in the last snapshot by numerically fitting to the (theoretical) solution 
( [2l| . |22]) with constant initial amplitude parameters vq and Av. The waves in the other 
snapshots were then shifted using only the theoretical velocity and the evolved time. The 
slope of the initial state was on purpose made much smaller than the correct one. The 
first three snapshots have been plotted with open circles, triangles and squares, and the 
remaining ones with solid circles. One can see that the initial state changes its slope and 
converges rapidly to the final common shape. It should be noted that this process does 
not generate any noticeable radiation, this is because the dissipation is strongest in the 
region which dominates the generation of radiation. 

The initial form ([H],[^) is actually very close to the final traveling wave solution. This 
is demonstrated in Fig. 3(b), where we have shown the difference between the theoretical 
approximation and the numerical result presented in Fig. 3(a). The deviation is only few 
percent of the original shape. 

If we change the sign of i n in ([22]), the SW starts to travel in the opposite direction as 
shown in Fig. 4. In this case, the profile of the wave is not constant, but rather its slope 
decreases gradually. 

The most interesting feature of the SW is how it adjusts itself when the initial velocity 
is different from that given by Eq. (§), this is illustrated in Fig. 5. If i n is smaller than 
the correct value (i.e., the velocity is too large, c.f. Eq. (|22|)), the system adjusts locally 
the upper level r + to agree with the given velocity (Fig. 5(a)), and this, in turn, gives rise 
to another SW at the top level. This new SW travels away from the original SW, but 
since this is the wrong direction this SW slowly flattens out. If, instead, we use initially 
z n (0) that is larger than the right one, the upper level of the wave increases (Fig. 5(b)), 
and the two SWs propagate in the opposite directions. For both the shocks, the velocities 
have the correct sign, therefore they quickly reach a stationary form without changing 
the slope (Fig. 5(b)). It is remarkable that all these major changes take place without 
any visible radiation. 
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4 Conclusions 



In this paper we have studied the existence of stationary SW's in the Toda lattice. Stable 
traveling waves of this type were found to exist when the friction in the lattice is of the 
"inner" type as per Eq. (|3[). This is not accidental. From Fig. 1 we see that, without 
dissipation, the system develops violent oscillations at the upper "shoulder" (level) of the 
SW. The second spatial derivative of the waveform has its maximum there as well, so that 
its time derivative provides a sufficiently strong balancing dissipative force according to 
Eq. (3). 

With the inner dissipation, the SW quickly attains its final stationary form. Using the 
continuum approximation, we have obtained the formula @ which relates the velocity of 
the SW to the asymptotic values of the lattice density in front of and beyond the SW. 
This formula proves to be very accurate even in the strongly discrete case. 

A very interesting feature of the SW is seen when the initial velocity and the asymptotic 
values do not obey Eq. The SW will then adjust itself to the velocity, and this always 
seems to happen through changing the upper level of the SW. This process gives rise to 
a new SW, which is traveling in the opposite direction. 

The existence of the stable SWs in the damped Toda lattice is a novel feature of the 
system and it is very likely that it may be observed in experiments. 
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Figure captions 



Fig. 1: (a) The time evolution of the initial state given by Eq. (|2l| , |22| ) in the absence of 
dissipation (7 = 0). The time interval is 100, vq = 1.0, Av = 1.0, k = 0.159. For 
better visibility we have left out the 9th and 10th snapshots. In this and all other 
figures we plot v n = exp(— r n ) — 1; (b) the same as Fig. 1(a) but with dissipation 
7 = 1.0. 

Fig. 2: The final waveforms of three numerical solutions with different values of the dissi- 
pation factor 7 (but with the same asymptotic densities obtained from v = 1.0, 
Av = 0.5). 

Fig. 3: (a) the approach to the stationary shape of the SW when the initial state has an 
incorrect slope (the time interval is 50, v Q = 1.0, Av = 3.0, k = 0.200, 7 = 0.5). 
Each curve has been shifted relative to the lattice in order to center them at the 
same lattice position (see the text). The three first curves are depicted, by open 
circles, triangles and squares, respectively, the remaining ones by solid circles; (b) 
the same as Fig. 3(a) but with the analytical approximation ( PT|) subtracted from 
each curve. 

Fig. 4: Snapshots of the lattice when the sign of initial i n (0) (hence the sign of V) is 
"incorrect", the initial state being in the rear (vq = 1.0, Av = 3.0, 7 = 0.5). 

Fig. 5: (a) The time evolution of an initial state with an "incorrect" value of the initial 
velocity (but with the correct sign) (t> = 1.0, Av = 3, 7 = 0.5). In this case, the 
"correct" velocity V is 1.809, while we set initially V = 10.0; (b): the same but 
with the initial velocity V = 1.0. 
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